Selective benefit of Glucosinolate and flavonoids

#Part 1 The purpose of this notebook is to determine the selective benefit of glucosinolate and flavonoid compounds through plant competition by creating selection gradients. Selection gradients will involve final body mass as the proxy for fitness, but this will be replaced with fitness once the measurement is in. Concentration will be on the x-axis. This analysis will account for family and greenhouse location.

#Part 2 The second purpose is to determine if glucosinolates and flavonoids influence suscpetibility to pathogens and if this susceptibility results in increased fitness

#Read in and prep data

library(ggplot2)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     

library(lme4)

Loading required package: Matrix

library(tidyr)

Attaching package: ‘tidyr’

The following objects are masked from ‘package:Matrix’:

    expand, pack, unpack

library(cowplot)

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************

library(grid)
library(gridExtra)
library(dplyr)

Attaching package: ‘dplyr’

The following object is masked from ‘package:gridExtra’:

    combine

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

source("GGPlot_Themes.R")

#read in data
rm(list=ls())
dat<-read.csv("DataSynthesis.csv")

#Remove maple controls data from the data set.
dat<-dat %>% filter(treatment!="mcnt") %>% droplevels()

#Assign family column
prefamily<-gsub("*.\\|","",dat$Tag)
dat$Family<-gsub("\\-.*","",prefamily)

#remove those with fertilizer treatment, and extra genotypes that are only in the alone treatment. 
dat<-dat[!grepl("i",dat$Sample,fixed=T),]

#Initializing columns to avoid constant error messages. 
dat$WhiteFungLogis<-NA
dat$BlackFungLogis<-NA

#Generating data frames with summarized leaf data (dat2) and summarized genotype data (dat3)

#Summarizing: Taking the mean value of leaves. This tibble contains data at the level of the plant means. 
dat2<-dat %>% group_by(Tag) %>% summarize(ChlorA=mean(ChlorA),ChlorB=mean(ChlorB),gluc_Conc=mean(gluc_Conc),flav_Conc=mean(flav_Conc),Family=first(Family),treatment=first(treatment),gh_row=first(gh_row),gh_bench=first(gh_bench),GM_TotalLeaf_Area=first(GM_TotalLeaf_Area),comp_number=first(comp_number),ThripsDam=mean(ThripsDam),WhiteFungDam=mean(WhiteFungDam),BlackPathDam=mean(BlackPathDam),Fern=mean(Fern),gh_col=first(gh_col))


#All of these genotypes died (15 Genotypes).(We are simply missing a final measurement for e|JBCHY1-1-50|Q|240) That is only 3% Mortality. 
dead<-dat2[is.na(dat2$GM_TotalLeaf_Area),]

#Removing those with dead competitors from the garlic mustard treatment. ("e|JBCHY1-1-50|Q|240") did not die, we are simply missing the final measurement for it.

dead_competitors<-dead %>% filter(treatment=="gm",Tag!="e|JBCHY1-1-50|Q|240") %>% select(comp_number)

#Removing those with dead competitors from the analysis. 
dat2<-dat2 %>% filter(!comp_number %in% dead_competitors$comp_number)


##Summarizing: Taking the mean value of plants. This tibble contains data at the level of the family means within each treatment. 
dat3<-dat2 %>% drop_na(GM_TotalLeaf_Area) %>% group_by(Family,treatment)  %>% summarize_if(is.numeric,mean)

#Searching for a fitness trade-off between the alone and interspecific competition treatment.

source("GGPlot_Themes.R")

#Calculating family means within treatment. 
TradeOff<-dat3 %>% filter(treatment=="a") %>% drop_na(GM_TotalLeaf_Area) %>%  select(LeafSizeAlone=GM_TotalLeaf_Area,Family,gluc_ConcAlone=gluc_Conc) %>% right_join(dat3 %>% filter(treatment=="m") %>%  select(LeafSizeMaple=GM_TotalLeaf_Area,Family,gluc_ConcMaple=gluc_Conc),by="Family")

#Calculating standard error of each family
StdErr<-dat2 %>% select(Family,treatment,GM_TotalLeaf_Area) %>% group_by(Family,treatment) %>% drop_na(GM_TotalLeaf_Area) %>%  summarize(StdErr=sd(GM_TotalLeaf_Area)/sqrt(length(GM_TotalLeaf_Area)),size=length(GM_TotalLeaf_Area))

#Shifting the data to be in long form. 
StdErr2<-StdErr %>% filter(treatment=="a") %>% select(StdErrAlone=StdErr,Family) %>% right_join(StdErr %>% filter(treatment=="m") %>% select(StdErrMaple=StdErr,Family),by="Family")

TradeOff2<-StdErr2 %>% left_join(TradeOff)

Joining, by = "Family"

#tiff("Selection_Figures/TradeOff.tiff", units="in", width=10, height=6, res=300)

ggplot(TradeOff2,aes(y=LeafSizeAlone,x=LeafSizeMaple))+
  geom_point()+
  geom_linerange(aes(ymin = LeafSizeAlone - StdErrAlone, 
                    ymax = LeafSizeAlone + StdErrAlone))+
  geom_errorbarh(aes(xmin = LeafSizeMaple - StdErrMaple,
                    xmax = LeafSizeMaple + StdErrMaple))+
theme_simple()+
  xlab("Performance with Maple")+
  ylab("Performance Alone")

#dev.off()
summary(lm(LeafSizeAlone~LeafSizeMaple,data=TradeOff2))

Call:
lm(formula = LeafSizeAlone ~ LeafSizeMaple, data = TradeOff2)

Residuals:
    Min      1Q  Median      3Q     Max 
-2019.5  -657.0  -265.5   907.7  2364.0 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9590.03571 1616.26626   5.933 6.87e-06 ***
LeafSizeMaple   -0.00936    0.18258  -0.051     0.96    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1244 on 21 degrees of freedom
Multiple R-squared:  0.0001251, Adjusted R-squared:  -0.04749 
F-statistic: 0.002628 on 1 and 21 DF,  p-value: 0.9596

#Visualizing genetic variation and greenhouse variation, which will be controlled for.

#GH Bench
ggplot(dat2)+
  geom_point(aes(y=gluc_Conc,x=gh_bench,colour=as.factor(gh_bench)))

#GH Col
ggplot(dat2)+
  geom_point(aes(y=gluc_Conc,x=gh_col,colour=as.factor(gh_bench)))

#Investigating genetic differences by treatment
#gluc_Conc
boxplot(gluc_Conc~Family,data=dat2[dat2$treatment=="a",])

boxplot(gluc_Conc~Family,data=dat2[dat2$treatment=="m",])

boxplot(gluc_Conc~Family,data=dat2[dat2$treatment=="gm",])

#bodymass
boxplot(GM_TotalLeaf_Area~Family,data=dat2[dat2$treatment=="a",])

boxplot(GM_TotalLeaf_Area~Family,data=dat2[dat2$treatment=="m",])

boxplot(GM_TotalLeaf_Area~Family,data=dat2[dat2$treatment=="gm",])

#Modelling: What influences performance (total leaf area)

#This is the full model.I expect that the influence of gluc conc and flav conc could vary between treatments because of allelopathy, therefore, i include interactions with these variables. However, i have no reason to think that pathogens or fern abundance could influence fitness differently in different treatments, therefore, no interactions are included.

#First lets ensure we are using the correct random effects.
fitfull0<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family), data=dat2)

fitfull<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench), data=dat2)

fitfull2<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench/gh_row), data=dat2)

fitfull3<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench/gh_col), data=dat2)

anova(fitfull0,fitfull) #Strong evidence to use bench.

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fitfull0: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
fitfull0:     WhiteFungDam + ThripsDam + Fern + (1 | Family)
fitfull: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
fitfull:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench)
         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
fitfull0   18 6603.2 6672.7 -3283.6   6567.2                        
fitfull    19 6576.9 6650.2 -3269.4   6538.9 28.37  1  1.002e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(fitfull,fitfull2) #Using bench/row is the same as the null model.

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fitfull: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
fitfull:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench)
fitfull2: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
fitfull2:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_row)
         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fitfull    19 6576.9 6650.2 -3269.4   6538.9                     
fitfull2   20 6577.7 6654.8 -3268.8   6537.7 1.2001  1     0.2733

anova(fitfull,fitfull3) #There is evidence to use bench/collumn however. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fitfull: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
fitfull:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench)
fitfull3: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
fitfull3:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col)
         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
fitfull    19 6576.9 6650.2 -3269.4   6538.9                      
fitfull3   20 6573.7 6650.9 -3266.9   6533.7 5.145  1    0.02331 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Modelling fixed effects. -----------

#Full Model
fitfull3<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench/gh_col), data=dat2)

#Removing three way interaction
fit.1<-update(fitfull3, ~.-treatment:gluc_Conc:flav_Conc)
anova(fitfull3,fit.1) #Good to remove

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.1: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.1:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.1:     treatment:gluc_Conc + treatment:flav_Conc + gluc_Conc:flav_Conc
fitfull3: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
fitfull3:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col)
         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.1      18 6573.6 6643.0 -3268.8   6537.6                     
fitfull3   20 6573.7 6650.9 -3266.9   6533.7 3.8462  2     0.1462

#Removing two way interaction gluc:flav
fit.2<-update(fit.1,~.-gluc_Conc:flav_Conc)
anova(fit.2,fit.1) #Good to remove. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.2: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.2:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.2:     treatment:gluc_Conc + treatment:flav_Conc
fit.1: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.1:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.1:     treatment:gluc_Conc + treatment:flav_Conc + gluc_Conc:flav_Conc
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.2   17 6573.6 6639.2 -3269.8   6539.6                     
fit.1   18 6573.6 6643.0 -3268.8   6537.6 2.0633  1     0.1509

#Removing treatment:flavonoid interactions
fit.3<-update(fit.2,~.-treatment:flav_Conc)
anova(fit.2,fit.3) #Good to remove. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.3: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.3:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.3:     treatment:gluc_Conc
fit.2: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.2:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.2:     treatment:gluc_Conc + treatment:flav_Conc
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.3   15 6570.3 6628.2 -3270.1   6540.3                     
fit.2   17 6573.6 6639.2 -3269.8   6539.6 0.6791  2     0.7121

fit.4<-update(fit.3,~.-treatment:gluc_Conc)
anova(fit.4,fit.3)  #That significantly affected the predictive power. Did Not Remove 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.4: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.4:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col)
fit.3: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.3:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.3:     treatment:gluc_Conc
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
fit.4   13 6572.4 6622.5 -3273.2   6546.4                       
fit.3   15 6570.3 6628.2 -3270.1   6540.3 6.0873  2    0.04766 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Flavonoid Concentration has missing samples. Therefore, I refit the model with and without flav_Conc, but using the same limited dataset.
#No flav_Conc
fit.4<-lmer(GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam +  
    ThripsDam + Fern +   treatment:gluc_Conc+ (1 | Family) + (1 | gh_bench/gh_col)  ,data=dat2[!is.na(dat2$flav_Conc),])
#Yes flav_Conc
fit.3<-lmer(GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam +  
    ThripsDam + Fern +   treatment:gluc_Conc+flav_Conc+ (1 | Family) + (1 | gh_bench/gh_col)  ,data=dat2[!is.na(dat2$flav_Conc),])

anova(fit.3,fit.4) #Flav_Conc is a significant predictor. 

refitting model(s) with ML (instead of REML)

Data: dat2[!is.na(dat2$flav_Conc), ]
Models:
fit.4: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.4:     ThripsDam + Fern + treatment:gluc_Conc + (1 | Family) + (1 | 
fit.4:     gh_bench/gh_col)
fit.3: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.3:     ThripsDam + Fern + treatment:gluc_Conc + flav_Conc + (1 | 
fit.3:     Family) + (1 | gh_bench/gh_col)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
fit.4   14 6574.0 6628.0 -3273.0   6546.0                       
fit.3   15 6570.3 6628.2 -3270.1   6540.3 5.7058  1    0.01691 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Now that I know flav_conc is significant, I will conduct the rest of the model simplification using the whole data set (i.e. flav_Conc excluded --fit4) 

fit.nf<-lmer(GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam +  
    ThripsDam + Fern +   treatment:gluc_Conc+ (1 | Family) + (1 | gh_bench/gh_col)  ,data=dat2)


fit.1.nf<-update(fit.nf,~.-BlackPathDam)
anova(fit.nf,fit.1.nf) #Black pathogen damage is a significant predictor. Did not remove.

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.1.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + WhiteFungDam + ThripsDam + 
fit.1.nf:     Fern + (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
fit.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.nf:     ThripsDam + Fern + treatment:gluc_Conc + (1 | Family) + (1 | 
fit.nf:     gh_bench/gh_col)
         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
fit.1.nf   13 6862.6 6913.3 -3418.3   6836.6                        
fit.nf     14 6854.1 6908.7 -3413.0   6826.1 10.487  1   0.001202 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.2.nf<-update(fit.nf,~.-WhiteFungDam)
anova(fit.nf,fit.2.nf)  #White Fungal Damage was not significant, however, both the AIC and BIC was lower with it included, suggesting it may be an important variable. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.2.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + ThripsDam + 
fit.2.nf:     Fern + (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
fit.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.nf:     ThripsDam + Fern + treatment:gluc_Conc + (1 | Family) + (1 | 
fit.nf:     gh_bench/gh_col)
         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.2.nf   13 6853.0 6903.7 -3413.5   6827.0                     
fit.nf     14 6854.1 6908.7 -3413.0   6826.1 0.8647  1     0.3524

dat2

fit.3.nf<-update(fit.nf,~.-ThripsDam)
anova(fit.nf,fit.3.nf)  #Thrips Damage was also not a significant predictor, however, both the AIC and BIC was lower with it included, suggesting it may be an important predictor. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.3.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.3.nf:     Fern + (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
fit.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.nf:     ThripsDam + Fern + treatment:gluc_Conc + (1 | Family) + (1 | 
fit.nf:     gh_bench/gh_col)
         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.3.nf   13 6853.6 6904.3 -3413.8   6827.6                     
fit.nf     14 6854.1 6908.7 -3413.0   6826.1 1.5362  1     0.2152

fit.4.nf<-update(fit.nf,~.-Fern)
anova(fit.nf,fit.4.nf)  #Fern abundance is significnat predictor performance.

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.4.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.4.nf:     ThripsDam + (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
fit.nf: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
fit.nf:     ThripsDam + Fern + treatment:gluc_Conc + (1 | Family) + (1 | 
fit.nf:     gh_bench/gh_col)
         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
fit.4.nf   13 6857.1 6907.8 -3415.5   6831.1                       
fit.nf     14 6854.1 6908.7 -3413.0   6826.1 4.9702  1    0.02579 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(lmerTest)

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

#The best model is therefore: 

#With limited dataset containing flavonoid data
summary(lmer(GM_TotalLeaf_Area ~ treatment*gluc_Conc+Fern+BlackPathDam +flav_Conc + (1 | gh_bench/gh_col) ,data=dat2))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: GM_TotalLeaf_Area ~ treatment * gluc_Conc + Fern + BlackPathDam +  
    flav_Conc + (1 | gh_bench/gh_col)
   Data: dat2

REML criterion at convergence: 6417.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2347 -0.5991  0.0668  0.5698  2.6996 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept)  548555   740.6  
 gh_bench        (Intercept) 1839408  1356.2  
 Residual                    7340803  2709.4  
Number of obs: 350, groups:  gh_col:gh_bench, 28; gh_bench, 5

Fixed effects:
                      Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           12088.20    2455.14   243.01   4.924 1.57e-06 ***
treatmentgm           -7338.20    2891.65   324.03  -2.538  0.01163 *  
treatmentm            -9033.07    3217.17   329.90  -2.808  0.00529 ** 
gluc_Conc             -4628.73    2410.11   323.99  -1.921  0.05567 .  
Fern                   -114.12      57.75   326.56  -1.976  0.04898 *  
BlackPathDam            -77.29      30.37   334.00  -2.545  0.01137 *  
flav_Conc              2630.10    1235.18   339.59   2.129  0.03395 *  
treatmentgm:gluc_Conc  3480.42    2920.34   322.79   1.192  0.23422    
treatmentm:gluc_Conc   8227.51    3344.05   330.58   2.460  0.01439 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm glc_Cn Fern   BlckPD flv_Cn trtmntg:_C
treatmentgm -0.770                                                       
treatmentm  -0.702  0.586                                                
gluc_Conc   -0.856  0.745   0.674                                        
Fern         0.012 -0.004  -0.068  -0.033                                
BlackPathDm -0.090  0.021   0.103  -0.015 -0.157                         
flav_Conc   -0.188  0.079   0.080  -0.265  0.032  0.139                  
trtmntgm:_C  0.754 -0.991  -0.575  -0.746 -0.002 -0.026 -0.058           
trtmntm:g_C  0.674 -0.562  -0.993  -0.654  0.060 -0.105 -0.079  0.558    

#with whole data set lacking flavonoid data
summary(lmer(GM_TotalLeaf_Area ~ treatment*gluc_Conc+Fern+BlackPathDam  + (1 | gh_bench/gh_col) ,data=dat2))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: GM_TotalLeaf_Area ~ treatment * gluc_Conc + Fern + BlackPathDam +  
    (1 | gh_bench/gh_col)
   Data: dat2

REML criterion at convergence: 6719.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2385 -0.6239  0.0654  0.5901  2.7601 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept)  366692   605.6  
 gh_bench        (Intercept) 1847886  1359.4  
 Residual                    7516712  2741.7  
Number of obs: 365, groups:  gh_col:gh_bench, 28; gh_bench, 5

Fixed effects:
                      Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           13372.80    2395.82   239.92   5.582 6.42e-08 ***
treatmentgm           -8127.48    2876.85   340.56  -2.825  0.00500 ** 
treatmentm            -8903.10    3119.91   346.11  -2.854  0.00458 ** 
gluc_Conc             -3558.56    2302.70   340.35  -1.545  0.12318    
Fern                   -132.91      57.26   326.48  -2.321  0.02089 *  
BlackPathDam            -86.08      30.01   352.48  -2.868  0.00438 ** 
treatmentgm:gluc_Conc  4235.27    2905.11   339.30   1.458  0.14580    
treatmentm:gluc_Conc   8010.78    3226.04   346.86   2.483  0.01349 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm glc_Cn Fern   BlckPD trtmntg:_C
treatmentgm -0.767                                                
treatmentm  -0.711  0.588                                         
gluc_Conc   -0.956  0.791   0.733                                 
Fern         0.014 -0.001  -0.048  -0.021                         
BlackPathDm -0.069  0.013   0.086   0.026 -0.152                  
trtmntgm:_C  0.752 -0.991  -0.578  -0.785 -0.006 -0.020           
trtmntm:g_C  0.683 -0.565  -0.993  -0.712  0.039 -0.087  0.562    

#The negative slope associated with the alone treamtment is no longer significant, suggesting that there is only a postive slope of glucosinolates in the maple treatment. However, the issue with this analysis is that is does not account for the competitive strength of the maple competitor, which may be correlated with glucosinolate concentration. 

#Model Diagnostics.

#GlucConcInteraction
fit.Whole<-lmer(GM_TotalLeaf_Area ~ treatment*gluc_Conc + BlackPathDam+Fern+flav_Conc+
(1 | Family) + (1 | gh_bench/gh_col) ,data=dat2)

confint(fit.Whole)

Computing profile confidence intervals ...
unexpected decrease in profile: using minstepnon-monotonic profile for (Intercept)bad spline fit for (Intercept): falling back to linear interpolation

                            2.5 %       97.5 %
.sig01                   274.5891  1300.474695
.sig02                     0.0000   893.955975
.sig03                   505.4262  2820.690738
.sigma                  2439.7371  2860.032525
(Intercept)             7363.1399 16799.855507
treatmentgm           -12701.8363 -1520.475152
treatmentm            -15041.5377 -2494.142955
gluc_Conc              -9383.2299   -53.492406
BlackPathDam            -143.4638   -24.193771
Fern                    -226.2257    -3.249272
flav_Conc                383.2981  5261.207752
treatmentgm:gluc_Conc  -2445.1744  8859.601582
treatmentm:gluc_Conc    1424.5156 14497.082917

coef(fit.Whole)

$`gh_col:gh_bench`
    (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam      Fern
1:1    12561.34   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
1:2    12098.72   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
1:3    11807.58   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
1:4    11587.20   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
2:1    12068.50   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
2:2    12218.84   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
2:3    12582.12   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
2:4    12628.25   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
2:5    12073.82   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
3:1    12744.74   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
3:2    12643.83   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
3:3    12768.44   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
3:4    12095.20   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
4:1    12091.12   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
4:2    12680.35   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
4:3    12499.17   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
4:4    11946.64   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
4:5    11806.81   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
5:1    11434.61   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
5:2    11600.05   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
5:3    11202.47   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
5:4    12179.72   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
5:5    12039.70   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
6:1    12007.78   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
6:2    10946.24   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
6:3    11566.32   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
6:4    11563.01   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
6:5    11922.11   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
    flav_Conc treatmentgm:gluc_Conc treatmentm:gluc_Conc
1:1  2832.853              3201.353             7947.925
1:2  2832.853              3201.353             7947.925
1:3  2832.853              3201.353             7947.925
1:4  2832.853              3201.353             7947.925
2:1  2832.853              3201.353             7947.925
2:2  2832.853              3201.353             7947.925
2:3  2832.853              3201.353             7947.925
2:4  2832.853              3201.353             7947.925
2:5  2832.853              3201.353             7947.925
3:1  2832.853              3201.353             7947.925
3:2  2832.853              3201.353             7947.925
3:3  2832.853              3201.353             7947.925
3:4  2832.853              3201.353             7947.925
4:1  2832.853              3201.353             7947.925
4:2  2832.853              3201.353             7947.925
4:3  2832.853              3201.353             7947.925
4:4  2832.853              3201.353             7947.925
4:5  2832.853              3201.353             7947.925
5:1  2832.853              3201.353             7947.925
5:2  2832.853              3201.353             7947.925
5:3  2832.853              3201.353             7947.925
5:4  2832.853              3201.353             7947.925
5:5  2832.853              3201.353             7947.925
6:1  2832.853              3201.353             7947.925
6:2  2832.853              3201.353             7947.925
6:3  2832.853              3201.353             7947.925
6:4  2832.853              3201.353             7947.925
6:5  2832.853              3201.353             7947.925

$Family
        (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam      Fern
BWPEM1     11892.98   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
CBMCK1     12266.13   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
CRSOSO     12294.71   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
CWRIC2     12114.65   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
DVGM       12208.20   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
EFCC2      11950.95   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
JARI1      12149.84   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
JBBLB2     11746.59   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
JBCHY1     11365.44   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
JWBOY      12441.57   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
KVEDG1     11739.19   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
MAVBEL2    12114.27   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
MHBUR1     12122.75   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
MHNAT1     11951.90   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
MSMID1     12144.45   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
PDVRT1     12163.71   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
RULEB1     12164.85   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
SMAKC1     12023.88   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
SMITH1     12551.35   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
VRCAN      11833.99   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
VRPET2     12018.79   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
VSGARO1    11938.44   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
WSSWM3     11922.34   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
        flav_Conc treatmentgm:gluc_Conc treatmentm:gluc_Conc
BWPEM1   2832.853              3201.353             7947.925
CBMCK1   2832.853              3201.353             7947.925
CRSOSO   2832.853              3201.353             7947.925
CWRIC2   2832.853              3201.353             7947.925
DVGM     2832.853              3201.353             7947.925
EFCC2    2832.853              3201.353             7947.925
JARI1    2832.853              3201.353             7947.925
JBBLB2   2832.853              3201.353             7947.925
JBCHY1   2832.853              3201.353             7947.925
JWBOY    2832.853              3201.353             7947.925
KVEDG1   2832.853              3201.353             7947.925
MAVBEL2  2832.853              3201.353             7947.925
MHBUR1   2832.853              3201.353             7947.925
MHNAT1   2832.853              3201.353             7947.925
MSMID1   2832.853              3201.353             7947.925
PDVRT1   2832.853              3201.353             7947.925
RULEB1   2832.853              3201.353             7947.925
SMAKC1   2832.853              3201.353             7947.925
SMITH1   2832.853              3201.353             7947.925
VRCAN    2832.853              3201.353             7947.925
VRPET2   2832.853              3201.353             7947.925
VSGARO1  2832.853              3201.353             7947.925
WSSWM3   2832.853              3201.353             7947.925

$gh_bench
  (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam      Fern
1    13929.35   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
2    11729.79   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
3    12457.09   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
4    11155.54   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
5    10971.92   -7099.007  -8758.159 -4720.332    -83.21952 -115.3353
  flav_Conc treatmentgm:gluc_Conc treatmentm:gluc_Conc
1  2832.853              3201.353             7947.925
2  2832.853              3201.353             7947.925
3  2832.853              3201.353             7947.925
4  2832.853              3201.353             7947.925
5  2832.853              3201.353             7947.925

attr(,"class")
[1] "coef.mer"

#no heteroscedasticity
plot(fit.Whole)

#fairly normal. 
qqnorm(resid(fit.Whole))

#The model appears to be a good fit.

#Visualization – The conditional benefit of glucosinolates (allelopathy)

aloneSlope<-function(x){
  y=-3966.60*x+ 13751.96
  return(y)
}
mapleSlope<-function(x){
  y=(-4254.2+8417.95)*x+13751.96 -9243.71
  return(y)
}

mustardSlope<-function(x){
  y=+13751.96-8376.67#Non significant slope
  return(y)
}

minM<-min(dat2$gluc_Conc[dat2$treatment=="m"],na.rm = T)
maxM<-max(dat2$gluc_Conc[dat2$treatment=="m"],na.rm = T)

minA<-min(dat2$gluc_Conc[dat2$treatment=="a"],na.rm = T)
maxA<-max(dat2$gluc_Conc[dat2$treatment=="a"],na.rm = T)

minG<-min(dat2$gluc_Conc[dat2$treatment=="gm"],na.rm = T)
maxG<-max(dat2$gluc_Conc[dat2$treatment=="gm"],na.rm = T)

#tiff("Selection_Figures/Gluc_Benefit.tiff", units="in", width=10, height=6, res=300)
library(ggplot2)
ggplot(dat2)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc,colour=treatment),size=2)+
  geom_segment(x=minA,xend=maxA,y=aloneSlope(minA),yend=aloneSlope(maxA),colour="#009E73",size=1.5)+
    geom_segment(x=minM,xend=maxM,y=mapleSlope(minM),yend=mapleSlope(maxM),colour="#E69F00",size=1.5)+
  geom_segment(x=minG,xend=maxG,y=mustardSlope(minG),yend=mustardSlope(maxG),colour="#56B4E9",size=1.5)+theme_simple()+
  ylab(bquote(bold("Performance\n(Total Leaf Area "~(mm^2)~")")))+xlab(bquote(bold("[Total Glucosinolate] " (mg/ml))))+
  scale_colour_manual(values=c("#009E73","#56B4E9","#E69F00"),labels=c("Alone","Garlic Mustard","Maple"))

#dev.off()

minA

[1] 0.7713421

#Visualizing the effect of flavonoid on fitness.

summary(lmer(GM_TotalLeaf_Area~treatment+flav_Conc+BlackPathDam+Fern+(1|Family)+(1|gh_bench/gh_col), data=dat2))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: GM_TotalLeaf_Area ~ treatment + flav_Conc + BlackPathDam + Fern +  
    (1 | Family) + (1 | gh_bench/gh_col)
   Data: dat2

REML criterion at convergence: 6473.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1990 -0.6278  0.0560  0.6004  2.6495 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept)  660632   812.8  
 Family          (Intercept)  210327   458.6  
 gh_bench        (Intercept) 1767108  1329.3  
 Residual                    7165525  2676.8  
Number of obs: 350, groups:  gh_col:gh_bench, 28; Family, 23; gh_bench, 5

Fixed effects:
             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)   7682.44    1156.80    33.30   6.641 1.42e-07 ***
treatmentgm  -3897.18     385.80   333.76 -10.102  < 2e-16 ***
treatmentm   -1133.35     364.02   330.79  -3.113  0.00201 ** 
flav_Conc     2471.19    1040.31   341.75   2.375  0.01808 *  
BlackPathDam   -77.69      30.19   335.79  -2.573  0.01050 *  
Fern          -125.75      57.56   330.33  -2.185  0.02962 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm flv_Cn BlckPD
treatmentgm -0.294                              
treatmentm  -0.260  0.475                       
flav_Conc   -0.795  0.198   0.137               
BlackPathDm -0.163 -0.043   0.021   0.090       
Fern        -0.035 -0.037  -0.073   0.030 -0.156

aloneSlope<-function(x){
  y=2024.62 *x+  8035.30
  return(y)
}
mapleSlope<-function(x){
  y=(2024.62 )*x+ 8035.30 -1067.02
  return(y)
}

mustardSlope<-function(x){
  y=2024.62 *x+ 8035.30-3746.91#Non significant slope
  return(y)
}

minM<-min(dat2$flav_Conc[dat2$treatment=="m"],na.rm = T)
maxM<-max(dat2$flav_Conc[dat2$treatment=="m"],na.rm = T)

minA<-min(dat2$flav_Conc[dat2$treatment=="a"],na.rm = T)
maxA<-max(dat2$flav_Conc[dat2$treatment=="a"],na.rm = T)

minG<-min(dat2$flav_Conc[dat2$treatment=="gm"],na.rm = T)
maxG<-max(dat2$flav_Conc[dat2$treatment=="gm"],na.rm = T)


#tiff("Selection_Figures/Flav_Benefit.tiff", units="in", width=10, height=6, res=300)
ggplot(dat2)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=flav_Conc,colour=treatment),size=2)+
  geom_segment(x=minA,xend=maxA,y=aloneSlope(minA),yend=aloneSlope(maxA),colour="#009E73",size=1.5)+
    geom_segment(x=minM,xend=maxM,y=mapleSlope(minM),yend=mapleSlope(maxM),colour="#E69F00",size=1.5)+theme_simple()+
  geom_segment(x=minG,xend=maxG,y=mustardSlope(minG),yend=mustardSlope(maxG),colour="#56B4E9",size=1.5)+theme_simple()+
  ylab(bquote(bold("Performance\n(Total Leaf Area "~(mm^2)~")")))+xlab(bquote(bold("[Total Flavonoid] " (mg/ml))))+
  scale_colour_manual(values=c("#009E73","#56B4E9","#E69F00"),labels=c("Alone","Garlic Mustard","Maple"))

#dev.off()

#Visualization– The Detriment of pathogens and ferns

summary(lmer(GM_TotalLeaf_Area ~ Fern+
(1 | Family) + (1 | gh_bench/gh_col) ,data=dat2))

boundary (singular) fit: see ?isSingular

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: GM_TotalLeaf_Area ~ Fern + (1 | Family) + (1 | gh_bench/gh_col)
   Data: dat2

REML criterion at convergence: 9607.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.85397 -0.70508  0.04291  0.74619  3.11787 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept)    96624  310.8  
 Family          (Intercept)        0    0.0  
 gh_bench        (Intercept)  2947085 1716.7  
 Residual                    10094168 3177.1  
Number of obs: 507, groups:  gh_col:gh_bench, 28; Family, 23; gh_bench, 5

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept) 7896.911    807.136    4.071   9.784 0.000561 ***
Fern        -226.020     56.097  432.082  -4.029 6.61e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
Fern -0.055
convergence code: 0
boundary (singular) fit: see ?isSingular

a<-ggplot(dat2)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=BlackPathDam))+theme_simple_multiCol()+
  geom_abline(intercept=8281.76,slope = -105.28,size=1.5)+
  xlab("Black Pathogen Damage")

b<-ggplot(dat2[dat2$WhiteFungDam<30,])+
  geom_point(aes(y=GM_TotalLeaf_Area,x=WhiteFungDam),colour="#999999")+theme_simple_multiCol()+
  theme(axis.title.x = element_text(color = "#999999", size = 16, face = "bold",margin=margin(3,0,3,0)),
      )+xlab("Powdery Mildew Damage")

c<-ggplot(dat2)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=ThripsDam),colour="#E69F00")+theme_simple_multiCol()+xlab("Thrips Damage")+
    theme(axis.title.x = element_text(color = "#E69F00", size = 16, face = "bold",margin=margin(3,0,3,0)))

d<-ggplot(dat2)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=Fern),colour="#009E73")+theme_simple_multiCol()+xlab("Fern Abundance")+
   geom_abline(intercept=8003.44,slope = -179.47,size=1.5,color="#009E73")+
  theme(axis.title.x = element_text(color = "#009E73", size = 16, face = "bold",margin=margin(3,0,3,0)))
   


plot<-plot_grid(a, b,ncol=2,rel_widths = c(1,1))

Removed 61 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).

plot2<-plot_grid(d, c,ncol=2,rel_widths = c(1,1))

Removed 9 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).

plot3<-plot_grid(a,b,c,d,ncol=2,rel_widths = c(1,1))

Removed 61 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).Removed 9 rows containing missing values (geom_point).

plot4<-plot_grid(a,b,c,ncol=1,rel_widths = c(1,1,1))

Removed 61 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).

plot5<-plot_grid(a,b,c,ncol=3,rel_widths = c(1,1,1))

Removed 61 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).Removed 61 rows containing missing values (geom_point).

y.grob<-textGrob(bquote(bold("Shoot Area "(mm^2))),gp=gpar(fontface="bold",fontsize=20),rot=90)


#tiff("Selection_Figures/PathogenEffect.tiff", units="in", width=14, height=6, res=300)
grid.arrange(plot,left=y.grob)

#dev.off()

#tiff("Selection_Figures/PathogenEffect2.tiff", units="in", width=14, height=6, res=300)
grid.arrange(plot2,left=y.grob)

#dev.off()

#tiff("Selection_Figures/PathogenEffect3.tiff", units="in", width=14, height=10, res=300)
grid.arrange(plot3,left=y.grob)

#dev.off

grid.arrange(plot4,left=y.grob)

#tiff("Selection_Figures/PathogenEffect5.tiff", units="in", width=16, height=5, res=300)
grid.arrange(plot5,left=y.grob)

#dev.off

#tiff("Selection_Figures/FernEffect.tiff", units="in", width=8, height=6, res=300)
grid.arrange(d,left=y.grob)

#dev.off

#—————————————– #Part 2

#Influence of gluc and flav on defence.

hist(dat$Fern)

hist(dat$ThripsDam)

hist(dat$WhiteFungDam)

hist(dat$BlackPathDam)

#This data is very zero inflated and a poisson model will be too overdispersed. I am afraid that a negative binomial will be as well. 

#Significant testing

#WhitePathDam – logistic regression

#Because of how zero inflated white pathogen damage is, i will use a logistic regression to model it.
dat2$WhiteFungLogis<-NA
dat2$WhiteFungLogis[dat2$WhiteFungDam==0]<-0
dat2$WhiteFungLogis[dat2$WhiteFungDam>0]<-1


#This is the biggest model that could converge. ... interaction with flavonoid could not.
fit_full_g<-glmer(WhiteFungLogis~treatment*gluc_Conc+flav_Conc+(1|Family),family=binomial,data=dat2[!is.na(dat2$flav_Conc),])

fit.1<-update(fit_full_g,~.-flav_Conc)
anova(fit_full_g,fit.1) #Flavonoid Concentration is not significant. 

Data: dat2[!is.na(dat2$flav_Conc), ]
Models:
fit.1: WhiteFungLogis ~ treatment + gluc_Conc + (1 | Family) + treatment:gluc_Conc
fit_full_g: WhiteFungLogis ~ treatment * gluc_Conc + flav_Conc + (1 | Family)
           npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.1         7 402.61 429.62 -194.31   388.61                     
fit_full_g    8 403.89 434.76 -193.95   387.89 0.7204  1      0.396

#Testing glucosinolate treatment interaction. 
fit.2<-glmer(WhiteFungLogis~treatment*gluc_Conc+(1|Family),family=binomial,data=dat2)

Model failed to converge with max|grad| = 0.00272048 (tol = 0.002, component 1)

fit.3<-glmer(WhiteFungLogis~treatment+gluc_Conc+(1|Family),family=binomial,data=dat2)
anova(fit.2,fit.3) #Glucosinolate:Treatment interaction is not significant.

Data: dat2
Models:
fit.3: WhiteFungLogis ~ treatment + gluc_Conc + (1 | Family)
fit.2: WhiteFungLogis ~ treatment * gluc_Conc + (1 | Family)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.3    5 416.98 436.48 -203.49   406.98                     
fit.2    7 416.63 443.93 -201.32   402.63 4.3439  2      0.114

#Testing glucosinolate involvment at all.
fit.4<-glmer(WhiteFungLogis~treatment+(1|Family),family=binomial,data=dat2[!is.na(dat2$gluc_Conc),])
fit.3<-glmer(WhiteFungLogis~treatment+gluc_Conc+(1|Family),family=binomial,data=dat2[!is.na(dat2$gluc_Conc),])
anova(fit.4,fit.3) #Glucosinolate Concentration is not a significant predictor

Data: dat2[!is.na(dat2$gluc_Conc), ]
Models:
fit.4: WhiteFungLogis ~ treatment + (1 | Family)
fit.3: WhiteFungLogis ~ treatment + gluc_Conc + (1 | Family)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.4    4 415.04 430.64 -203.52   407.04                     
fit.3    5 416.98 436.48 -203.49   406.98 0.0644  1     0.7997

#Testing effect of treatment
fit.4<-glmer(WhiteFungLogis~treatment+(1|Family),family=binomial,data=dat2)
fit.5<-glmer(WhiteFungLogis~1+(1|Family),family=binomial,data=dat2)
anova(fit.4,fit.5)

Data: dat2
Models:
fit.5: WhiteFungLogis ~ 1 + (1 | Family)
fit.4: WhiteFungLogis ~ treatment + (1 | Family)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
fit.5    2 523.93 532.17 -259.97   519.93                       
fit.4    4 519.23 535.71 -255.62   511.23 8.7022  2    0.01289 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#treatment is significant......

summary(fit.4) #Garlic mustard in the maple treatment have less occurence of fungal colonization. 

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: WhiteFungLogis ~ treatment + (1 | Family)
   Data: dat2

     AIC      BIC   logLik deviance df.resid 
   519.2    535.7   -255.6    511.2      451 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3885 -0.6239 -0.4587  0.9947  3.3329 

Random effects:
 Groups Name        Variance Std.Dev.
 Family (Intercept) 0.544    0.7376  
Number of obs: 455, groups:  Family, 23

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.6835     0.2395  -2.854  0.00431 **
treatmentgm  -0.4101     0.2613  -1.569  0.11656   
treatmentm   -0.8037     0.2756  -2.916  0.00355 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg
treatmentgm -0.509        
treatmentm  -0.474  0.447 

plot(fit.4)

#Visualizing – effect of treatment on proportion of fungal abundance.

#Summarizing for Display: generating frequency of white funal infection by treatment
plot<-dat2 %>% drop_na(WhiteFungLogis) %>% group_by(treatment) %>% summarize(PercWhitFung=sum(WhiteFungLogis)/length(WhiteFungLogis)*100)


#tiff("Defence_Figures/TreatMeanWhiteFung.tiff", units="in", width=8, height=5, res=300)
ggplot(plot)+
  geom_col(aes(x=treatment,y=PercWhitFung,fill=treatment))+theme_simple()+ylab("Fungal Abundance\n (% Infected)")+
  scale_y_continuous(breaks = seq(0,30,5))+
scale_x_discrete(name="",labels=c("Alone","Garlic Mustard","Maple"))+
  scale_fill_manual(values=c("#009E73","#56B4E9","#E69F00"),labels=c("Alone","Garlic Mustard","Maple"))+
  theme_simple_multiCol()+theme(axis.title.y =  element_text(color = "black", size = 16, face = "bold",margin=margin(3,20,3,0)))

#dev.off()

#Visualizing — the distribution of pathogens by treatment.

plot2<-dat2 %>% drop_na(BlackPathDam) %>% group_by(treatment) %>% summarize(BlackPathAve=mean(BlackPathDam,na.rm=T))

plot3<-dat2 %>% drop_na(ThripsDam) %>% group_by(treatment) %>% summarize(ThripsDam=mean(ThripsDam,na.rm=T))

plot4<-dat2 %>% drop_na(Fern) %>% group_by(treatment) %>% summarize(Fern=mean(Fern,na.rm=T))


#Black Pathogen abundance by treamtent
ggplot(plot2)+
  geom_col(aes(x=treatment,y=BlackPathAve,fill=treatment))+theme_simple()+ylab("Fungal Abundance\n (% colonized)")+xlab("Treatment")+theme(legend.position = "none")

#Thrips abundance by treatment. 
ggplot(plot3)+
  geom_col(aes(x=treatment,y=ThripsDam,fill=treatment))+theme_simple()+ylab("Fungal Abundance\n (% colonized)")+xlab("Treatment")+theme(legend.position = "none")

#I would need to account for the garlic mustard being psuedoreplicated iwth the fern data, because each gm in the gm treatment was grown with another gm in the same pot. 
ggplot(plot4)+
  geom_col(aes(x=treatment,y=Fern,fill=treatment))+theme_simple()+ylab("Average Fern Abundance (#Ferns/pot)")+xlab("Treatment")+theme(legend.position = "none")+
  scale_x_discrete(name="",labels=c("Alone","Garlic Mustard","Maple"))+
  scale_fill_manual(values=c("#009E73","#56B4E9","#E69F00"),labels=c("Alone","Garlic Mustard","Maple"))+
  theme_simple_multiCol()+theme(axis.title.y =  element_text(color = "black", size = 16, face = "bold",margin=margin(3,20,3,0)))

#dev.off()

#Modelling – what influences White pathogen damage?

#Rounding to integers for white fung dam. 
dat$WhiteFungDam<-ceiling(dat$WhiteFungDam)

fit_full<-glmer(WhiteFungDam~treatment*gluc_Conc+treatment*flav_Conc+(1|Family/Tag)+(1|gh_bench),family=poisson,data=dat)

Model failed to converge with max|grad| = 1.93052 (tol = 0.002, component 1)

fit.1<-update(fit_full,~.-flav_Conc:treatment)

Model failed to converge with max|grad| = 0.368802 (tol = 0.002, component 1)

anova(fit.1,fit_full) #Fit1 is a better model

Data: dat
Models:
fit.1: WhiteFungDam ~ treatment + gluc_Conc + flav_Conc + (1 | Family/Tag) + 
fit.1:     (1 | gh_bench) + treatment:gluc_Conc
fit_full: WhiteFungDam ~ treatment * gluc_Conc + treatment * flav_Conc + 
fit_full:     (1 | Family/Tag) + (1 | gh_bench)
         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
fit.1      10 1011.1 1056.3 -495.53   991.06                    
fit_full   12 1016.3 1070.6 -496.15   992.29     0  2          1

fit.2<-update(fit.1,~.-gluc_Conc:treatment)

boundary (singular) fit: see ?isSingular

anova(fit.2,fit.1) #models are the same... eliminating interaction.

Data: dat
Models:
fit.2: WhiteFungDam ~ treatment + gluc_Conc + flav_Conc + (1 | Family/Tag) + 
fit.2:     (1 | gh_bench)
fit.1: WhiteFungDam ~ treatment + gluc_Conc + flav_Conc + (1 | Family/Tag) + 
fit.1:     (1 | gh_bench) + treatment:gluc_Conc
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.2    8 1010.9 1047.1 -497.43   994.85                     
fit.1   10 1011.1 1056.3 -495.53   991.06 3.7917  2     0.1502

fit.3<-update(fit.2,~.-gluc_Conc)

boundary (singular) fit: see ?isSingular

anova(fit.3,fit.2) #The model including glucConc is better. 

Data: dat
Models:
fit.3: WhiteFungDam ~ treatment + flav_Conc + (1 | Family/Tag) + (1 | 
fit.3:     gh_bench)
fit.2: WhiteFungDam ~ treatment + gluc_Conc + flav_Conc + (1 | Family/Tag) + 
fit.2:     (1 | gh_bench)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
fit.3    7 1019.2 1050.9 -502.60  1005.21                        
fit.2    8 1010.9 1047.1 -497.43   994.85 10.356  1   0.001291 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.2<-update(fit.2,data=dat[!is.na(dat$flav_Conc),]) #Using a model with a data subset.

boundary (singular) fit: see ?isSingular

fit.3<-update(fit.2,~.-flav_Conc) 

boundary (singular) fit: see ?isSingular

anova(fit.3,fit.2)#flav_Conc is unimportant, even though the AIC and BIC are lower.

Data: dat[!is.na(dat$flav_Conc), ]
Models:
fit.3: WhiteFungDam ~ treatment + gluc_Conc + (1 | Family/Tag) + (1 | 
fit.3:     gh_bench)
fit.2: WhiteFungDam ~ treatment + gluc_Conc + flav_Conc + (1 | Family/Tag) + 
fit.2:     (1 | gh_bench)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
fit.3    7 1009.4 1041.1 -497.71   995.42                     
fit.2    8 1010.9 1047.1 -497.43   994.85 0.5708  1       0.45

#Therefore, the best model is:
summary(glmer(WhiteFungDam~treatment+gluc_Conc+(1|Family/Tag)+(1|gh_bench),family=poisson,data=dat))

boundary (singular) fit: see ?isSingular

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: WhiteFungDam ~ treatment + gluc_Conc + (1 | Family/Tag) + (1 |  
    gh_bench)
   Data: dat

     AIC      BIC   logLik deviance df.resid 
  1033.7   1065.6   -509.9   1019.7      697 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0689 -0.2817 -0.2273 -0.1792  2.4492 

Random effects:
 Groups     Name        Variance Std.Dev.
 Tag:Family (Intercept) 3.5607   1.8870  
 Family     (Intercept) 0.1535   0.3918  
 gh_bench   (Intercept) 0.0000   0.0000  
Number of obs: 704, groups:  Tag:Family, 433; Family, 23; gh_bench, 5

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.19132    0.74202  -0.258  0.79653   
treatmentgm  0.03378    0.33151   0.102  0.91883   
treatmentm  -0.81258    0.35172  -2.310  0.02087 * 
gluc_Conc   -2.32106    0.71153  -3.262  0.00111 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm
treatmentgm -0.306                
treatmentm  -0.341  0.484         
gluc_Conc   -0.907  0.126   0.175 
convergence code: 0
boundary (singular) fit: see ?isSingular

#the inclusion of GH_Bench does not affect the outcome of the model, thereofre, I am excluding it as the other model was saturated and took alot time to compute. 
summary(glmer(WhiteFungDam~treatment+gluc_Conc+(1|Family/Tag),family=poisson,data=dat))

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: WhiteFungDam ~ treatment + gluc_Conc + (1 | Family/Tag)
   Data: dat

     AIC      BIC   logLik deviance df.resid 
  1031.7   1059.1   -509.9   1019.7      698 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0689 -0.2816 -0.2273 -0.1792  2.4492 

Random effects:
 Groups     Name        Variance Std.Dev.
 Tag:Family (Intercept) 3.5609   1.8870  
 Family     (Intercept) 0.1535   0.3918  
Number of obs: 704, groups:  Tag:Family, 433; Family, 23

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.19135    0.74206  -0.258  0.79651   
treatmentgm  0.03379    0.33152   0.102  0.91882   
treatmentm  -0.81257    0.35173  -2.310  0.02088 * 
gluc_Conc   -2.32108    0.71157  -3.262  0.00111 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm
treatmentgm -0.306                
treatmentm  -0.341  0.484         
gluc_Conc   -0.907  0.126   0.175 

#Using a Permutation test to determine the probability of getting a significant glucosinolate concentration predictor, if the data is ramdomly sampled. 

datPerTest<-dat
pStore<-c()
for(i in 1:100){
  #Randomize glucosinolate concentration.
  datPerTest$gluc_Conc<-sample(dat$gluc_Conc,length(dat$gluc_Conc),replace = F)
  
  #Run Model with randomized glucConc and extract p-value
  glucPval<-summary(glmer(WhiteFungDam~treatment+gluc_Conc+(1|Family/Tag),family=poisson,data=datPerTest))$coef[3,4]
  
  #Store P value.
  pStore[i]<-glucPval
}

boundary (singular) fit: see ?isSingular
unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 2 negative eigenvaluesvariance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RXvariance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RXboundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular
Model failed to converge with max|grad| = 0.115676 (tol = 0.002, component 1)Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?boundary (singular) fit: see ?isSingular
Model failed to converge with max|grad| = 0.11074 (tol = 0.002, component 1)Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?boundary (singular) fit: see ?isSingular
boundary (singular) fit: see ?isSingular

#Visuallizing – white fungal damage is influenced by gluc and treatmnet.

PoisSlope=function(x,int){
  y=exp(-2.32292*x+int)
  return(y)
}

glucplot=seq(min(dat$gluc_Conc),max(dat$gluc_Conc),length.out = 695)

FungA<-PoisSlope(glucplot,-0.15701)
FungM<-PoisSlope(glucplot,-0.15701-0.82060)
FungGM<-PoisSlope(glucplot,-0.15701+0.02483)

fit.g<-glm(WhiteFungDam~gluc_Conc+treatment,family=poisson,data=dat)
summary(fit.g)

predict(fit.g,type="response",newdata=data.frame("treatment"=rep("a",695),"gluc_Conc"=glucplot))

datv<-dat[!is.na(dat$WhiteFungDam),]

#tiff("Defence_Figures/GlucosinolateFungi.tiff", units="in", width=10, height=6, res=300)
ggplot(datv[datv$WhiteFungDam<30,])+
  geom_point(aes(y=WhiteFungDam,x=gluc_Conc,colour=treatment))+
    scale_colour_manual(values=c("#009E73","#56B4E9","#E69F00","#009E73","#56B4E9","#E69F00"),labels=c("Alone","Garlic Mustard","Maple"))+theme_simple()+
  scale_y_continuous(breaks=c(0,2,4,6,8,10,12))+
  geom_path(x=glucplot,y=FungA,colour="#009E73")+
  geom_path(x=glucplot,y=FungM,colour="#E69F00")+
  ylab("Powdery Mildew Damage")+
xlab(bquote(bold("[Total Glucosinolate] " (mg/ml))))
#dev.off()

#Modelling – What influences Black Pathogen Damage?

sum(pStore<0.05)/length(pStore)

[1] 0.56

#Visuallizing – Black Pathogen damage is influenced by flavonols.

PoisSlope=function(x){
  y=exp(-1.1150*x+1.6816)
  return(y)
}


flavplot=seq(min(dat$flav_Conc,na.rm = T),max(dat$flav_Conc,na.rm=T),length.out = 687)

flavy<-PoisSlope(flavplot)

#tiff("Defence_Figures/FlavonoidBlackPath.tiff", units="in", width=10, height=6, res=300)
ggplot(dat[!is.na(dat$flav_Conc)&!is.na(dat$flav_Conc),])+
  geom_point(aes(y=BlackPathDam,x=flav_Conc))+theme_simple()+
  geom_path(x=flavplot,y=flavy,size=1,colour="#999999")+
  
  scale_y_continuous(breaks=c(0,5,10,15,20,25,30))+
  ylab("Black Pathogen Damage")+
xlab(bquote(bold("[Total Flavonoid] " (mg/ml))))
#dev.off()

#What influences Thrips damage?

sum(pStore<0.05)/length(pStore)

[1] 0.47

#Visualizing - influence of flavonoids on thrips damage.

PoisSlope=function(x){
  y=exp(-2.3319*x+3.5249)
  return(y)
}


flavplot=seq(min(dat$flav_Conc,na.rm = T),max(dat$flav_Conc,na.rm=T),length.out = 687)

flavy<-PoisSlope(flavplot)

#tiff("Defence_Figures/FlavonoidThrips.tiff", units="in", width=10, height=6, res=300)
ggplot(dat[!is.na(dat$flav_Conc)&!is.na(dat$flav_Conc),])+
  geom_point(aes(y=ThripsDam,x=flav_Conc))+theme_simple()+
  geom_path(x=flavplot,y=flavy,size=1,colour="#999999")+
  
  scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40))+
  ylab("Thrips Damage")+
xlab(bquote(bold("[Total Flavonoid] " (mg/ml))))
#dev.off()

#Visualizing – effect of flavonoids on fern abundance.

PoisSlope=function(x,int){
  y=exp((-0.02353 -3.12475)*x+int)
  return(y)
}
exp(-0.6571)

flavplot=seq(min(dat$flav_Conc,na.rm = T),max(dat$flav_Conc,na.rm=T),length.out = 707)

flavyA<-PoisSlope(flavplot,-2.1683)
flavyM<-PoisSlope(flavplot,-1.60546-2.81450)
flavyGM<-PoisSlope(flavplot,-2.1683-0.2786)


#tiff("Defence_Figures/FlavonoidFern.tiff", units="in", width=10, height=6, res=300)
ggplot(dat)+
  geom_point(aes(y=Fern,x=flav_Conc,colour=treatment))+theme_simple()+
 # geom_path(x=flavplot,y=flavyA,size=1,colour="#009E73")
  #geom_path(x=flavplot,y=flavyGM,size=1,colour="#56B4E9")
  geom_path(x=flavplot,y=flavyM,size=1,colour="#E69F00")+
      scale_colour_manual(values=c("#009E73","#56B4E9","#E69F00"),labels=c("Alone","Garlic Mustard","Maple"))+
  scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40))+
  ylab("Fern Abundance")+
xlab(bquote(bold("[Total Flavonoid] " (mg/ml))))
#dev.off()

How can we know that healthy plants dont just exhibit more secondary compounds and not that those with more secondary compounds are healthier?